Usenet 1993 July

home *** CD-ROM | disk | FTP | other *** search

/ Usenet 1993 July / InfoMagic USENET CD-ROM July 1993.ISO / sources / misc / volume11 / ephem4.12 / part03 < prev next >

Wrap

Text File | 1990-03-10 | 59.6 KB | 2,161 lines

Newsgroups: comp.sources.misc subject: v11i004: ephem, 3 of 7 From: ecd@cs.umn.edu@ncs-med.UUCP (Elwood C. Downey) Sender: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc) Posting-number: Volume 11, Issue 4 Submitted-by: ecd@cs.umn.edu@ncs-med.UUCP (Elwood C. Downey) Archive-name: ephem4.12/part03 # This is the first line of a "shell archive" file. # This means it contains several files that can be extracted into # the current directory when run with the sh shell, as follows: # sh < this_file_name # This is file 3. echo x mainmenu.c sed -e 's/^X//' << 'EOFxEOF' > mainmenu.c X/* printing routines for the main (upper) screen. X */ X X#include <stdio.h> X#include <math.h> X#include "astro.h" X#include "circum.h" X#include "screen.h" X X/* #define PC_GRAPHICS */ X#ifdef PC_GRAPHICS X#define JOINT 207 X#define VERT 179 X#define HORIZ 205 X#else X#define JOINT '-' X#define VERT '|' X#define HORIZ '-' X#endif X Xmm_borders() X{ X char line[NC+1], *lp; X register i; X X lp = line; X for (i = 0; i < NC; i++) X *lp++ = HORIZ; X *lp = '\0'; X f_string (R_PLANTAB-1, 1, line); X for (i = R_TOP; i < R_PLANTAB-1; i++) X f_char (i, COL2-2, VERT); X f_char (R_PLANTAB-1, COL2-2, JOINT); X for (i = R_TOP; i < R_PLANTAB-1; i++) X f_char (i, COL3-2, VERT); X f_char (R_PLANTAB-1, COL3-2, JOINT); X for (i = R_LST; i < R_PLANTAB-1; i++) X f_char (i, COL4-2, VERT); X f_char (R_PLANTAB-1, COL4-2, JOINT); X} X X/* print the permanent labels, on the top menu and the planet names X * in the bottom section. X */ Xmm_labels() X{ X f_string (R_TZN, C_TZN, "LT"); X f_string (R_UT, C_UT, "UTC"); X f_string (R_JD, C_JD, "JulianDate"); X f_string (R_WATCH, C_WATCH, "Watch"); X f_string (R_SRCH, C_SRCH, "Search"); X f_string (R_PLOT, C_PLOT, "Plot"); X f_string (R_ALTM, C_ALTM, "Menu"); X X f_string (R_LST, C_LST, "LST"); X f_string (R_DAWN, C_DAWN, "Dawn"); X f_string (R_DUSK, C_DUSK, "Dusk"); X f_string (R_LON, C_LON, "NiteLn"); X f_string (R_NSTEP, C_NSTEP, "NStep"); X f_string (R_STPSZ, C_STPSZ, "StpSz"); X X f_string (R_LAT, C_LAT, "Lat"); X f_string (R_LONG, C_LONG, "Long"); X f_string (R_HEIGHT, C_HEIGHT, "Elev"); X f_string (R_TEMP, C_TEMP, "Temp"); X f_string (R_PRES, C_PRES, "AtmPr"); X f_string (R_TZONE, C_TZONE, "TZ"); X f_string (R_EPOCH, C_EPOCH, "Epoch"); X X f_string (R_PLANTAB, C_OBJ, "Ob"); X f_string (R_SUN, C_OBJ, "Su"); X f_string (R_MOON, C_OBJ, "Mo"); X f_string (R_MERCURY, C_OBJ, "Me"); X f_string (R_VENUS, C_OBJ, "Ve"); X f_string (R_MARS, C_OBJ, "Ma"); X f_string (R_JUPITER, C_OBJ, "Ju"); X f_string (R_SATURN, C_OBJ, "Sa"); X f_string (R_URANUS, C_OBJ, "Ur"); X f_string (R_NEPTUNE, C_OBJ, "Ne"); X f_string (R_PLUTO, C_OBJ, "Pl"); X f_string (R_OBJX, C_OBJ, "X"); X} X X/* print all the time/date/where related stuff: the Now structure. X * print in a nice order, based on the field locations, as much as possible. X */ Xmm_now (np, all) XNow *np; Xint all; X{ X char buf[32]; X double lmjd = mjd - tz/24.0; X double jd = mjd + 2415020L; X double tmp; X X sprintf (buf, "%-3.3s", tznm); X f_string (R_TZN, C_TZN, buf); X f_time (R_LT, C_LT, mjd_hr(lmjd)); X f_date (R_LD, C_LD, mjd_day(lmjd)); X X f_time (R_UT, C_UTV, mjd_hr(mjd)); X f_date (R_UD, C_UD, mjd_day(mjd)); X X sprintf (buf, "%14.5f", jd); X (void) flog_log (R_JD, C_JDV, jd); X f_string (R_JD, C_JDV, buf); X X now_lst (np, &tmp); X f_time (R_LST, C_LSTV, tmp); X X if (all) { X f_gangle (R_LAT, C_LATV, lat); X f_gangle (R_LONG, C_LONGV, -lng); /* + west */ X X tmp = height * 2.093e7; /* want to see ft, not earth radii */ X sprintf (buf, "%5g ft", tmp); X (void) flog_log (R_HEIGHT, C_HEIGHTV, tmp); X f_string (R_HEIGHT, C_HEIGHTV, buf); X X tmp = 9./5.*temp + 32.0; /* want to see degrees F, not C */ X (void) flog_log (R_TEMP, C_TEMPV, tmp); X sprintf (buf, "%6g F", tmp); X f_string (R_TEMP, C_TEMPV, buf); X X tmp = pressure / 33.86; /* want to see in. Hg, not mBar */ X (void) flog_log (R_PRES, C_PRESV, tmp); X sprintf (buf, "%5.2f in", tmp); X f_string (R_PRES, C_PRESV, buf); X X f_signtime (R_TZONE, C_TZONEV, tz); X X if (epoch == EOD) X f_string (R_EPOCH, C_EPOCHV, "(OfDate)"); X else { X mjd_year (epoch, &tmp); X f_double (R_EPOCH, C_EPOCHV, "%8.1f", tmp); X } X } X X /* print the calendar for local day, if new month/year. */ X mm_calendar (np, all > 1); X} X X/* display dawn/dusk/length-of-night times. X */ Xmm_twilight (np, force) XNow *np; Xint force; X{ X double dusk, dawn; X double tmp; X int status; X X if (!twilight_cir (np, &dawn, &dusk, &status) && !force) X return; X X if (status != 0) { X f_blanks (R_DAWN, C_DAWNV, 5); X f_blanks (R_DUSK, C_DUSKV, 5); X f_string (R_LON, C_LONV, "-----"); X return; X } X X f_mtime (R_DAWN, C_DAWNV, dawn); X f_mtime (R_DUSK, C_DUSKV, dusk); X tmp = dawn - dusk; range (&tmp, 24.0); X f_mtime (R_LON, C_LONV, tmp); X} X Xmm_newcir (y) Xint y; X{ X static char ncmsg[] = "NEW CIRCUMSTANCES"; X static char nomsg[] = " "; X static int last_y = -1; X X if (y != last_y) { X f_string (R_NEWCIR, C_NEWCIR, y ? ncmsg : nomsg); X last_y = y; X } X} X Xstatic Xmm_calendar (np, force) XNow *np; Xint force; X{ X static char *mnames[] = { X "January", "February", "March", "April", "May", "June", X "July", "August", "September", "October", "November", "December" X }; X static int last_m, last_y; X static double last_tz; X char str[64]; X int m, y; X double d; X int f, nd; X int r; X double jd0; X X /* get local m/d/y. do nothing if still same month and not forced. */ X mjd_cal (mjd_day(mjd-tz/24.0), &m, &d, &y); X if (m == last_m && y == last_y && tz == last_tz && !force) X return; X last_m = m; X last_y = y; X last_tz = tz; X X /* find day of week of first day of month */ X cal_mjd (m, 1.0, y, &jd0); X mjd_dow (jd0, &f); X if (f < 0) { X /* can't figure it out - too hard before Gregorian */ X int i; X for (i = 8; --i >= 0; ) X f_string (R_CAL+i, C_CAL, " "); X return; X } X X /* print header */ X f_blanks (R_CAL, C_CAL, 20); X sprintf (str, "%s %4d", mnames[m-1], y); X f_string (R_CAL, C_CAL + (20 - (strlen(mnames[m-1]) + 5))/2, str); X f_string (R_CAL+1, C_CAL, "Su Mo Tu We Th Fr Sa"); X X /* find number of days in this month */ X mjd_dpm (jd0, &nd); X X /* print the calendar */ X for (r = 0; r < 6; r++) { X char row[7*3+1], *rp = row; X int c; X for (c = 0; c < 7; c++) { X int i = r*7+c; X if (i < f || i >= f + nd) X sprintf (rp, " "); X else X sprintf (rp, "%2d ", i-f+1); X rp += 3; X } X row[sizeof(row)-2] = '\0'; /* don't print last blank; causes wrap*/ X f_string (R_CAL+2+r, C_CAL, row); X } X X /* over print the new and full moons for this month. X * TODO: don't really know which dates to use here (see moonnf()) X * so try several to be fairly safe. have to go back to 4/29/1988 X * to find the full moon on 5/1 for example. X */ X mm_nfmoon (jd0-3, tz, m, f); X mm_nfmoon (jd0+15, tz, m, f); X} X Xstatic Xmm_nfmoon (jd, tzone, m, f) Xdouble jd, tzone; Xint m, f; X{ X static char nm[] = "NM", fm[] = "FM"; X double dm; X int mm, ym; X double jdn, jdf; X int di; X X moonnf (jd, &jdn, &jdf); X mjd_cal (jdn-tzone/24.0, &mm, &dm, &ym); X if (m == mm) { X di = dm + f - 1; X f_string (R_CAL+2+di/7, C_CAL+3*(di%7), nm); X } X mjd_cal (jdf-tzone/24.0, &mm, &dm, &ym); X if (m == mm) { X di = dm + f - 1; X f_string (R_CAL+2+di/7, C_CAL+3*(di%7), fm); X } X} EOFxEOF len=`wc -c < mainmenu.c` if expr $len != 6662 > /dev/null then echo Length of mainmenu.c is $len but it should be 6662. fi echo x moon.c sed -e 's/^X//' << 'EOFxEOF' > moon.c X#include <stdio.h> X#include <math.h> X#include "astro.h" X X/* given the mjd, find the geocentric ecliptic longitude, lam, and latitude, X * bet, and horizontal parallax, hp for the moon. X * N.B. series for long and lat are good to about 10 and 3 arcseconds. however, X * math errors cause up to 100 and 30 arcseconds error, even if use double. X * why?? suspect highly sensitive nature of difference used to get m1..6. X * N.B. still need to correct for nutation. then for topocentric location X * further correct for parallax and refraction. X */ Xmoon (mjd, lam, bet, hp) Xdouble mjd; Xdouble *lam, *bet, *hp; X{ X double t, t2; X double ld; X double ms; X double md; X double de; X double f; X double n; X double a, sa, sn, b, sb, c, sc, e, e2, l, g, w1, w2; X double m1, m2, m3, m4, m5, m6; X X t = mjd/36525.; X t2 = t*t; X X m1 = mjd/27.32158213; X m1 = 360.0*(m1-(long)m1); X m2 = mjd/365.2596407; X m2 = 360.0*(m2-(long)m2); X m3 = mjd/27.55455094; X m3 = 360.0*(m3-(long)m3); X m4 = mjd/29.53058868; X m4 = 360.0*(m4-(long)m4); X m5 = mjd/27.21222039; X m5 = 360.0*(m5-(long)m5); X m6 = mjd/6798.363307; X m6 = 360.0*(m6-(long)m6); X X ld = 270.434164+m1-(.001133-.0000019*t)*t2; X ms = 358.475833+m2-(.00015+.0000033*t)*t2; X md = 296.104608+m3+(.009192+.0000144*t)*t2; X de = 350.737486+m4-(.001436-.0000019*t)*t2; X f = 11.250889+m5-(.003211+.0000003*t)*t2; X n = 259.183275-m6+(.002078+.000022*t)*t2; X X a = degrad(51.2+20.2*t); X sa = sin(a); X sn = sin(degrad(n)); X b = 346.56+(132.87-.0091731*t)*t; X sb = .003964*sin(degrad(b)); X c = degrad(n+275.05-2.3*t); X sc = sin(c); X ld = ld+.000233*sa+sb+.001964*sn; X ms = ms-.001778*sa; X md = md+.000817*sa+sb+.002541*sn; X f = f+sb-.024691*sn-.004328*sc; X de = de+.002011*sa+sb+.001964*sn; X e = 1-(.002495+7.52e-06*t)*t; X e2 = e*e; X X ld = degrad(ld); X ms = degrad(ms); X n = degrad(n); X de = degrad(de); X f = degrad(f); X md = degrad(md); X X l = 6.28875*sin(md)+1.27402*sin(2*de-md)+.658309*sin(2*de)+ X .213616*sin(2*md)-e*.185596*sin(ms)-.114336*sin(2*f)+ X .058793*sin(2*(de-md))+.057212*e*sin(2*de-ms-md)+ X .05332*sin(2*de+md)+.045874*e*sin(2*de-ms)+.041024*e*sin(md-ms); X l = l-.034718*sin(de)-e*.030465*sin(ms+md)+.015326*sin(2*(de-f))- X .012528*sin(2*f+md)-.01098*sin(2*f-md)+.010674*sin(4*de-md)+ X .010034*sin(3*md)+.008548*sin(4*de-2*md)-e*.00791*sin(ms-md+2*de)- X e*.006783*sin(2*de+ms); X l = l+.005162*sin(md-de)+e*.005*sin(ms+de)+.003862*sin(4*de)+ X e*.004049*sin(md-ms+2*de)+.003996*sin(2*(md+de))+ X .003665*sin(2*de-3*md)+e*.002695*sin(2*md-ms)+ X .002602*sin(md-2*(f+de))+e*.002396*sin(2*(de-md)-ms)- X .002349*sin(md+de); X l = l+e2*.002249*sin(2*(de-ms))-e*.002125*sin(2*md+ms)- X e2*.002079*sin(2*ms)+e2*.002059*sin(2*(de-ms)-md)- X .001773*sin(md+2*(de-f))-.001595*sin(2*(f+de))+ X e*.00122*sin(4*de-ms-md)-.00111*sin(2*(md+f))+.000892*sin(md-3*de); X l = l-e*.000811*sin(ms+md+2*de)+e*.000761*sin(4*de-ms-2*md)+ X e2*.000704*sin(md-2*(ms+de))+e*.000693*sin(ms-2*(md-de))+ X e*.000598*sin(2*(de-f)-ms)+.00055*sin(md+4*de)+.000538*sin(4*md)+ X e*.000521*sin(4*de-ms)+.000486*sin(2*md-de); X l = l+e2*.000717*sin(md-2*ms); X *lam = ld+degrad(l); X range (lam, 2*PI); X X g = 5.12819*sin(f)+.280606*sin(md+f)+.277693*sin(md-f)+ X .173238*sin(2*de-f)+.055413*sin(2*de+f-md)+.046272*sin(2*de-f-md)+ X .032573*sin(2*de+f)+.017198*sin(2*md+f)+.009267*sin(2*de+md-f)+ X .008823*sin(2*md-f)+e*.008247*sin(2*de-ms-f); X g = g+.004323*sin(2*(de-md)-f)+.0042*sin(2*de+f+md)+ X e*.003372*sin(f-ms-2*de)+e*.002472*sin(2*de+f-ms-md)+ X e*.002222*sin(2*de+f-ms)+e*.002072*sin(2*de-f-ms-md)+ X e*.001877*sin(f-ms+md)+.001828*sin(4*de-f-md)-e*.001803*sin(f+ms)- X .00175*sin(3*f); X g = g+e*.00157*sin(md-ms-f)-.001487*sin(f+de)-e*.001481*sin(f+ms+md)+ X e*.001417*sin(f-ms-md)+e*.00135*sin(f-ms)+.00133*sin(f-de)+ X .001106*sin(f+3*md)+.00102*sin(4*de-f)+.000833*sin(f+4*de-md)+ X .000781*sin(md-3*f)+.00067*sin(f+4*de-2*md); X g = g+.000606*sin(2*de-3*f)+.000597*sin(2*(de+md)-f)+ X e*.000492*sin(2*de+md-ms-f)+.00045*sin(2*(md-de)-f)+ X .000439*sin(3*md-f)+.000423*sin(f+2*(de+md))+ X .000422*sin(2*de-f-3*md)-e*.000367*sin(ms+f+2*de-md)- X e*.000353*sin(ms+f+2*de)+.000331*sin(f+4*de); X g = g+e*.000317*sin(2*de+f-ms+md)+e2*.000306*sin(2*(de-ms)-f)- X .000283*sin(md+3*f); X w1 = .0004664*cos(n); X w2 = .0000754*cos(c); X *bet = degrad(g)*(1-w1-w2); X X *hp = .950724+.051818*cos(md)+.009531*cos(2*de-md)+.007843*cos(2*de)+ X .002824*cos(2*md)+.000857*cos(2*de+md)+e*.000533*cos(2*de-ms)+ X e*.000401*cos(2*de-md-ms)+e*.00032*cos(md-ms)-.000271*cos(de)- X e*.000264*cos(ms+md)-.000198*cos(2*f-md); X *hp = *hp+.000173*cos(3*md)+.000167*cos(4*de-md)-e*.000111*cos(ms)+ X .000103*cos(4*de-2*md)-.000084*cos(2*md-2*de)- X e*.000083*cos(2*de+ms)+.000079*cos(2*de+2*md)+.000072*cos(4*de)+ X e*.000064*cos(2*de-ms+md)-e*.000063*cos(2*de+ms-md)+ X e*.000041*cos(ms+de); X *hp = *hp+e*.000035*cos(2*md-ms)-.000033*cos(3*md-2*de)- X .00003*cos(md+de)-.000029*cos(2*(f-de))-e*.000029*cos(2*md+ms)+ X e2*.000026*cos(2*(de-ms))-.000023*cos(2*(f-de)+md)+ X e*.000019*cos(4*de-ms-md); X *hp = degrad(*hp); X} EOFxEOF len=`wc -c < moon.c` if expr $len != 5143 > /dev/null then echo Length of moon.c is $len but it should be 5143. fi echo x moonnf.c sed -e 's/^X//' << 'EOFxEOF' > moonnf.c X#include <stdio.h> X#include <math.h> X#include "astro.h" X X#define unw(w,z) ((w)-floor((w)/(z))*(z)) X X/* given a modified Julian date, mjd, return the mjd of the new X * and full moons about then, mjdn and mjdf. X * TODO: exactly which ones does it find? eg: X * 5/28/1988 yields 5/15 and 5/31 X * 5/29 6/14 6/29 X */ Xmoonnf (mjd, mjdn, mjdf) Xdouble mjd; Xdouble *mjdn, *mjdf; X{ X int mo, yr; X double dy; X double mjd0; X double k, tn, tf, t; X X mjd_cal (mjd, &mo, &dy, &yr); X cal_mjd (1, 0., yr, &mjd0); X k = (yr-1900+((mjd-mjd0)/365))*12.3685; X k = floor(k+0.5); X tn = k/1236.85; X tf = (k+0.5)/1236.85; X t = tn; X m (t, k, mjdn); X t = tf; X k += 0.5; X m (t, k, mjdf); X} X Xstatic Xm (t, k, mjd) Xdouble t, k; Xdouble *mjd; X{ X double t2, a, a1, b, b1, c, ms, mm, f, ddjd; X X t2 = t*t; X a = 29.53*k; X c = degrad(166.56+(132.87-9.173e-3*t)*t); X b = 5.8868e-4*k+(1.178e-4-1.55e-7*t)*t2+3.3e-4*sin(c)+7.5933E-1; X ms = 359.2242+360*unw(k/1.236886e1,1)-(3.33e-5+3.47e-6*t)*t2; X mm = 306.0253+360*unw(k/9.330851e-1,1)+(1.07306e-2+1.236e-5*t)*t2; X f = 21.2964+360*unw(k/9.214926e-1,1)-(1.6528e-3+2.39e-6*t)*t2; X ms = unw(ms,360); X mm = unw(mm,360); X f = unw(f,360); X ms = degrad(ms); X mm = degrad(mm); X f = degrad(f); X ddjd = (1.734e-1-3.93e-4*t)*sin(ms)+2.1e-3*sin(2*ms) X -4.068e-1*sin(mm)+1.61e-2*sin(2*mm)-4e-4*sin(3*mm) X +1.04e-2*sin(2*f)-5.1e-3*sin(ms+mm)-7.4e-3*sin(ms-mm) X +4e-4*sin(2*f+ms)-4e-4*sin(2*f-ms)-6e-4*sin(2*f+mm) X +1e-3*sin(2*f-mm)+5e-4*sin(ms+2*mm); X a1 = (long)a; X b = b+ddjd+(a-a1); X b1 = (long)b; X a = a1+b1; X b = b-b1; X *mjd = a + b; X} EOFxEOF len=`wc -c < moonnf.c` if expr $len != 1557 > /dev/null then echo Length of moonnf.c is $len but it should be 1557. fi echo x nutation.c sed -e 's/^X//' << 'EOFxEOF' > nutation.c X#include <stdio.h> X#include <math.h> X#include "astro.h" X X/* given the modified JD, mjd, find the nutation in obliquity, *deps, and X * the nutation in longitude, *dpsi, each in radians. X */ Xnutation (mjd, deps, dpsi) Xdouble mjd; Xdouble *deps, *dpsi; X{ X static double lastmjd, lastdeps, lastdpsi; X double ls, ld; /* sun's mean longitude, moon's mean longitude */ X double ms, md; /* sun's mean anomaly, moon's mean anomaly */ X double nm; /* longitude of moon's ascending node */ X double t, t2; /* number of Julian centuries of 36525 days since X * Jan 0.5 1900. X */ X double tls, tnm, tld; /* twice above */ X double a, b; /* temps */ X X if (mjd == lastmjd) { X *deps = lastdeps; X *dpsi = lastdpsi; X return; X } X X t = mjd/36525.; X t2 = t*t; X X a = 100.0021358*t; X b = 360.*(a-(long)a); X ls = 279.697+.000303*t2+b; X X a = 1336.855231*t; X b = 360.*(a-(long)a); X ld = 270.434-.001133*t2+b; X X a = 99.99736056000026*t; X b = 360.*(a-(long)a); X ms = 358.476-.00015*t2+b; X X a = 13255523.59*t; X b = 360.*(a-(long)a); X md = 296.105+.009192*t2+b; X X a = 5.372616667*t; X b = 360.*(a-(long)a); X nm = 259.183+.002078*t2-b; X X /* convert to radian forms for use with trig functions. X */ X tls = 2*degrad(ls); X nm = degrad(nm); X tnm = 2*degrad(nm); X ms = degrad(ms); X tld = 2*degrad(ld); X md = degrad(md); X X /* find delta psi and eps, in arcseconds. X */ X lastdpsi = (-17.2327-.01737*t)*sin(nm)+(-1.2729-.00013*t)*sin(tls) X +.2088*sin(tnm)-.2037*sin(tld)+(.1261-.00031*t)*sin(ms) X +.0675*sin(md)-(.0497-.00012*t)*sin(tls+ms) X -.0342*sin(tld-nm)-.0261*sin(tld+md)+.0214*sin(tls-ms) X -.0149*sin(tls-tld+md)+.0124*sin(tls-nm)+.0114*sin(tld-md); X lastdeps = (9.21+.00091*t)*cos(nm)+(.5522-.00029*t)*cos(tls) X -.0904*cos(tnm)+.0884*cos(tld)+.0216*cos(tls+ms) X +.0183*cos(tld-nm)+.0113*cos(tld+md)-.0093*cos(tls-ms) X -.0066*cos(tls-nm); X X /* convert to radians. X */ X lastdpsi = degrad(lastdpsi/3600); X lastdeps = degrad(lastdeps/3600); X X lastmjd = mjd; X *deps = lastdeps; X *dpsi = lastdpsi; X} EOFxEOF len=`wc -c < nutation.c` if expr $len != 2010 > /dev/null then echo Length of nutation.c is $len but it should be 2010. fi echo x objx.c sed -e 's/^X//' << 'EOFxEOF' > objx.c X/* functions to save the user-definable object. X * this way, once defined, the object can be quieried for position just like the X * other bodies, with objx_cir(). X */ X X#include <math.h> X#include "astro.h" X#include "circum.h" X#include "screen.h" X Xstatic int onflag; /* !=0 while object x is active */ X Xstatic int xtype; /* see flags, below */ X#define FIXED 0 /* just simple ra/dec object */ X#define ELLIPTICAL 1 /* elliptical orbital elements */ X#define PARABOLIC 2 /* parabolic " */ X X/* the defining info about object x. X */ Xtypedef struct { X double f_ra; /* ra, rads, at given epoch */ X double f_dec; /* dec, rads, at given epoch */ X double f_epoch; /* the given epoch, as an mjd */ X} ObjF; /* fixed object */ Xtypedef struct { X double e_ep; /* epoch of perihelion, as an mjd */ X double e_p; /* orbital period, in (tropical?) years */ X double e_pl; /* perihelion longitude, degrees */ X double e_q; /* eccentricity */ X double e_inc; /* inclination, degrees */ X double e_om; /* longitude of ascending node, degrees */ X double e_sma; /* semi-major axis, AU */ X double e_m1, e_m2; /* magnitude model coefficients */ X} ObjE; /* object in heliocentric elliptical orbit */ Xtypedef struct { X double p_ep; /* epoch of perihelion, as an mjd */ X double p_inc; /* inclination, rads */ X double p_qp; /* perihelion distance, AU */ X double p_ap; /* argument of perihelion, rads. */ X double p_om; /* longitude of ascending node, rads */ X double p_epoch; /* reference epoch, as an mjd */ X double p_m1, p_m2; /* magnitude model coefficients */ X} ObjP; /* object in heliocentric parabolic trajectory */ X Xstatic ObjF objf; Xstatic ObjE obje; Xstatic ObjP objp; X X/* run when Objx is picked from menu. X * let op define object and turn it on and off. X */ Xobjx_setup() X{ X static char *p[4] = { X "Fixed coordinates", "Elliptical elements", "Parabolic elements" X }; X int f; X X /* start pointing at selection for current object type */ X switch (xtype) { X case FIXED: f = 0; break; X case ELLIPTICAL: f = 1; break; X case PARABOLIC: f = 2; break; X } X X ask: X p[3] = onflag ? "On" : "Off"; X switch (f = popup (p, f, 4)) { X case 0: objx_dfixed((char**)0); goto ask; X case 1: objx_delliptical((char**)0); goto ask; X case 2: objx_dparabolic((char**)0); goto ask; X case 3: onflag ^= 1; break; X } X} X X/* turn "on" or "off" but don't forget facts about object x. X */ Xobjx_on () X{ X onflag = 1; X} X Xobjx_off() X{ X onflag = 0; X} X X/* return true if object is now on, else 0. X */ Xobjx_ison() X{ X return (onflag); X} X X/* fill in info about object x. X * most arguments and conditions are the same as for plans(). X * only difference is that mag is already apparent, not absolute magnitude. X * this is called by body_cir() for object x just like plans() is called X * for the planets. X */ Xobjx_cir (jd, lpd0, psi0, rp0, rho0, lam, bet, mag) Xdouble jd; /* mjd now */ Xdouble *lpd0; /* heliocentric longitude, or NOHELIO */ Xdouble *psi0; /* heliocentric latitude, or 0 if *lpd0 set to NOHELIO */ Xdouble *rp0; /* distance from the sun, or 0 */ Xdouble *rho0; /* true distance from the Earth, or 0 */ Xdouble *lam; /* apparent geocentric ecliptic longitude */ Xdouble *bet; /* apparent geocentric ecliptic latitude */ Xdouble *mag; /* APPARENT magnitude, or NOMAG */ X{ X switch (xtype) { X case FIXED: { X double xr, xd; X xr = objf.f_ra; X xd = objf.f_dec; X if (objf.f_epoch != jd) X precess (objf.f_epoch, jd, &xr, &xd); X eq_ecl (jd, xr, xd, bet, lam); X X *lpd0 = NOHELIO; X *psi0 = *rp0 = *rho0 = 0.0; X *mag = NOMAG; X } X break; X case PARABOLIC: { X double inc, ap, om; X double lpd, psi, rp, rho; X double dt; X int pass; X /* two passes to correct lam and bet for light travel time. */ X dt = 0.0; X for (pass = 0; pass < 2; pass++) { X reduce_elements (objp.p_epoch, jd-dt, objp.p_inc, objp.p_ap, X objp.p_om, &inc, &ap, &om); X comet (jd-dt, objp.p_ep, inc, ap, objp.p_qp, om, X &lpd, &psi, &rp, &rho, lam, bet); X if (pass == 0) { X *lpd0 = lpd; X *psi0 = psi; X *rp0 = rp; X *rho0 = rho; X } X dt = rho*5.775518e-3; /* au to light-days */ X } X *mag = objp.p_m1 + 5*log10(*rho0) + objp.p_m2*log10(*rp0); X } X break; X case ELLIPTICAL: { X /* this is basically the same code as pelement() and plans() X * combined and simplified for the special case of osculating X * (unperturbed) elements. X */ X double a0, a1; X double t, aa, ml, dayl, dt, lg, lsn, rsn; X double nu, ea; X double ma, rp, lp, om, lo, slo, clo; X double inc, psi, spsi, cpsi; X double y, lpd, rpd, ll, rho, sll, cll; X int pass; X X a0 = obje.e_pl - 360.0/36525.0*100.0*obje.e_ep/obje.e_p; X a1 = 100.0/obje.e_p; X X t = jd/36525.0; X aa = a1*t; X ml = a0 + 360.0*(aa - (long)aa); X range (&ml, 360.0); X dayl = a1*9.856263e-3; X X dt = 0; X sunpos (jd, &lsn, &rsn); X lg = lsn + PI; X X for (pass = 0; pass < 2; pass++) { X ma = degrad(ml-obje.e_pl-dt*dayl); X anomaly (ma, obje.e_q, &nu, &ea); X rp = obje.e_sma*(1-obje.e_q*obje.e_q)/(1+obje.e_q*cos(nu)); X lp = raddeg(nu)+obje.e_pl; X lp = degrad(lp); X om = degrad(obje.e_om); X lo = lp-om; X slo = sin(lo); X clo = cos(lo); X inc = degrad(obje.e_inc); X spsi = slo*sin(inc); X y = slo*cos(inc); X psi = asin(spsi); X lpd = atan(y/clo)+om; X if (clo<0) lpd += PI; X range (&lpd, 2*PI); X cpsi = cos(psi); X rpd = rp*cpsi; X ll = lpd-lg; X rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll)); X dt = rho*5.775518e-3; X if (pass == 0) { X *lpd0 = lpd; X *psi0 = psi; X *rp0 = rp; X *rho0 = rho; X } X } X X sll = sin(ll); X cll = cos(ll); X *lam = atan(rsn*sll/(rpd-rsn*cll))+lpd; X /* *lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI; */ X range (lam, 2*PI); X *bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*rsn*sll)); X *mag = obje.e_m1 + 5*log10(*rho0) + obje.e_m2*log10(*rp0); X } X break; X } X} X X/* define objx based on the ephem.cfg line, s. X * the "OBJX " keyword is already skipped. X * format: xtype,[other fields, as per corresponding ObjX typedef] X * N.B. we replace all ',' with '\0' within s IN PLACE. X * return 0 if ok, else print reason why not with f_msg() and return -1. X */ Xobjx_define (s) Xchar *s; X{ X#define MAXARGS 20 X char *av[MAXARGS]; /* point to each field for easy reference */ X char c; X int ac; X X /* parse into comma separated fields */ X ac = 0; X av[0] = s; X do { X c = *s++; X if (c == ',' || c == '\0') { X s[-1] = '\0'; X av[++ac] = s; X } X } while (c); X X if (ac < 1) { X f_msg ("No type for OBJX"); X return (-1); X } X X switch (av[0][0]) { X case 'f': X if (ac != 4) { X f_msg ("Need ra,dec,epoch for \"fixed\" OBJX"); X return (-1); X } X objx_dfixed (av+1); X break; X case 'p': X if (ac != 9) { X f_msg ("Need ep,inc,ap,qp,om,epoch,abs,lum for \"parabolic\" OBJX"); X return (-1); X } X objx_dparabolic (av+1); X break; X case 'e': X if (ac != 10) { X f_msg ("Need ep,p,pl,ecc,inc,om,sma,abs,lum for \"elliptical\" OBJX"); X return (-1); X } X objx_delliptical (av+1); X break; X default: X f_msg ("Unknown OBJX type"); X return (-1); X } X X return (0); X} X X/* define a fixed object. X * args in av, in order, are ra, dec, end reference epoch. X * if av then it is a list of strings to use for each parameter, else must X * ask for each. the av option is for cracking the ephem.cfg line. X * if asking show current settings and leave unchanged if hit RETURN. X * END aborts without making any more changes. X * xtype is set to FIXED. X * N.B. we don't error check av in any way, not even for length. X */ Xstatic Xobjx_dfixed (av) Xchar *av[]; X{ X char buf[NC]; X char *bp; X int sts; X X xtype = FIXED; X X if (av) { X bp = av[0]; X sts = 1; X } else { X static char p[] = "RA (h:m:s): ("; X f_prompt (p); X f_ra (R_PROMPT, C_PROMPT+sizeof(p)-1, objf.f_ra); X printf (") "); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) { X int h, m, s; X f_dec_sexsign (radhr(objf.f_ra), &h, &m, &s); X f_sscansex (bp, &h, &m, &s); X sex_dec (h, m, s, &objf.f_ra); X objf.f_ra = hrrad(objf.f_ra); X } X X if (av) { X bp = av[1]; X sts = 1; X } else { X static char p[] = "Dec (d:m:s): ("; X f_prompt (p); X f_angle (R_PROMPT, C_PROMPT+sizeof(p)-1, objf.f_dec); X printf (") "); X sts = read_line (buf, 9+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) { X int dg, m, s; X f_dec_sexsign (raddeg(objf.f_dec), &dg, &m, &s); X f_sscansex (bp, &dg, &m, &s); X sex_dec (dg, m, s, &objf.f_dec); X objf.f_dec = degrad(objf.f_dec); X } X X if (av) { X bp = av[2]; X sts = 1; X } else { X static char p[] = "Reference epoch (UT Date, m/d.d/y or year.d): "; X double y; X f_prompt (p); X mjd_year (objf.f_epoch, &y); X printf ("(%g) ", y); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X crack_year (bp, &objf.f_epoch); X} X X/* define an object in an elliptical heliocentric orbit. X * args in av, in order, are epoch of perihelion, orbital period, perihelion X * longitude, eccentricity, inclination, longitude of ascending node, and X * semi-major axis. X * if av then it is a list of strings to use for each parameter, else must X * ask for each. the av option is for cracking the ephem.cfg line. X * if asking show current settings and leave unchanged if hit RETURN. X * END aborts without making any more changes. X * xtype is set to ELLIPTICAL. X * N.B. we don't error check av in any way, not even for length. X */ Xstatic Xobjx_delliptical(av) Xchar *av[]; X{ X char buf[NC]; X char *bp; X int sts; X X xtype = ELLIPTICAL; X X if (av) { X bp = av[0]; X sts = 1; X } else { X static char p[] = X "Epoch of perihelion (UT Date, m/d.d/y or year.d): "; X int m, y; X double d; X f_prompt (p); X mjd_cal (obje.e_ep, &m, &d, &y); X printf ("(%d/%g/%d) ", m, d, y); X sts = read_line (buf, 15); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X crack_year (bp, &obje.e_ep); X X if (av) { X bp = av[1]; X sts = 1; X } else { X static char p[] = "Oribital period (years):"; X f_prompt (p); X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_p); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X obje.e_p = atof(bp); X X if (av) { X bp = av[2]; X sts = 1; X } else { X static char p[] = "Perihelion longitude (degs):"; X f_prompt (p); X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_pl); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X obje.e_pl = atof(bp); X X if (av) { X bp = av[3]; X sts = 1; X } else { X static char p[] = "Eccentricity:"; X f_prompt (p); X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_q); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X obje.e_q = atof(bp); X X if (av) { X bp = av[4]; X sts = 1; X } else { X static char p[] = "Inclination (degs):"; X f_prompt (p); X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_inc); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X obje.e_inc = atof(bp); X X if (av) { X bp = av[5]; X sts = 1; X } else { X static char p[] = "Longitude of ascending node (degs):"; X f_prompt (p); X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_om); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X obje.e_om = atof(bp); X X if (av) { X bp = av[6]; X sts = 1; X } else { X static char p[] = "Semi-major axis (AU):"; X f_prompt (p); X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_sma); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X obje.e_sma = atof(bp); X X if (av) { X bp = av[7]; X sts = 1; X } else { X static char p[] = "Absolute magnitude:"; X f_prompt (p); X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_m1); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X obje.e_m1 = atof(bp); X X if (av) { X bp = av[8]; X sts = 1; X } else { X static char p[] = "Luminosity index:"; X f_prompt (p); X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_m2); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X obje.e_m2 = atof(bp); X} X X/* define an object in heliocentric parabolic orbit. X * args in av, in order, are epoch of perihelion, inclination, argument of X * perihelion, perihelion distance, longitude of ascending node and X * reference epoch. X * if av then it is a list of strings to use for each parameter, else must X * ask for each. the av option is for cracking the ephem.cfg line. X * if asking show current settings and leave unchanged if hit RETURN. X * END aborts without making any more changes. X * xtype is set to PARABOLIC. X * N.B. we don't error check av in any way, not even for length. X */ Xstatic Xobjx_dparabolic(av) Xchar *av[]; X{ X char buf[NC]; X char *bp; X int sts; X X xtype = PARABOLIC; X X if (av) { X bp = av[0]; X sts = 1; X } else { X static char p[] = X "Epoch of perihelion (UT Date, m/d.d/y or year.d): "; X int m, y; X double d; X f_prompt (p); X mjd_cal (objp.p_ep, &m, &d, &y); X printf ("(%d/%g/%d) ", m, d, y); X sts = read_line (buf, PW-sizeof(p)); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X crack_year (bp, &objp.p_ep); X X if (av) { X bp = av[1]; X sts = 1; X } else { X static char p[] = "Inclination (degs):"; X f_prompt (p); X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ",raddeg(objp.p_inc)); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X objp.p_inc = degrad(atof(bp)); X X if (av) { X bp = av[2]; X sts = 1; X } else { X static char p[] = "Argument of perihelion (degs):"; X f_prompt (p); X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", raddeg(objp.p_ap)); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X objp.p_ap = degrad(atof(bp)); X X if (av) { X bp = av[3]; X sts = 1; X } else { X static char p[] = "Perihelion distance (AU):"; X f_prompt (p); X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", objp.p_qp); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X objp.p_qp = atof(bp); X X if (av) { X bp = av[4]; X sts = 1; X } else { X static char p[] = "Longitude of ascending node (degs):"; X f_prompt (p); X f_double(R_PROMPT,C_PROMPT+sizeof(p), "(%g) ", raddeg(objp.p_om)); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X objp.p_om = degrad(atof(bp)); X X if (av) { X bp = av[5]; X sts = 1; X } else { X static char p[] = "Reference epoch (UT Date, m/d.d/y or year.d): "; X double y; X f_prompt (p); X mjd_year (objp.p_epoch, &y); X printf ("(%g) ", y); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X crack_year (bp, &objp.p_epoch); X X if (av) { X bp = av[6]; X sts = 1; X } else { X static char p[] = "Absolute magnitude:"; X f_prompt (p); X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", objp.p_m1); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X objp.p_m1 = atof(bp); X X if (av) { X bp = av[7]; X sts = 1; X } else { X static char p[] = "Luminosity index:"; X f_prompt (p); X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", objp.p_m2); X sts = read_line (buf, 8+1); X if (sts < 0) X return; X bp = buf; X } X if (sts > 0) X objp.p_m2 = atof(bp); X} X X/* given either a decimal year (xxxx. something) or a calendar (x/x/x) X * convert it to an mjd and store it at *p; X */ Xstatic Xcrack_year (bp, p) Xchar *bp; Xdouble *p; X{ X if (decimal_year(bp)) { X double y = atof (bp); X year_mjd (y, p); X } else { X int m, y; X double d; X mjd_cal (objp.p_ep, &m, &d, &y); /* init with current */ X f_sscandate (bp, &m, &d, &y); X cal_mjd (m, d, y, p); X } X} EOFxEOF len=`wc -c < objx.c` if expr $len != 16243 > /dev/null then echo Length of objx.c is $len but it should be 16243. fi echo x obliq.c sed -e 's/^X//' << 'EOFxEOF' > obliq.c X#include <stdio.h> X#include "astro.h" X X/* given the modified Julian date, mjd, find the obliquity of the X * ecliptic, *eps, in radians. X */ Xobliquity (mjd, eps) Xdouble mjd; Xdouble *eps; X{ X static double lastmjd, lasteps; X X if (mjd != lastmjd) { X double t; X t = mjd/36525.; X lasteps = degrad(2.345229444E1 - ((((-1.81E-3*t)+5.9E-3)*t+4.6845E1)*t)/3600.0); X lastmjd = mjd; X } X *eps = lasteps; X} EOFxEOF len=`wc -c < obliq.c` if expr $len != 409 > /dev/null then echo Length of obliq.c is $len but it should be 409. fi echo x parallax.c sed -e 's/^X//' << 'EOFxEOF' > parallax.c X#include <stdio.h> X#include <math.h> X#include "astro.h" X X/* given true ha and dec, tha and tdec, the geographical latitude, phi, the X * height above sea-level (as a fraction of the earths radius, 6378.16km), X * ht, and the equatorial horizontal parallax, ehp, find the apparent X * ha and dec, aha and adec allowing for parallax. X * all angles in radians. ehp is the angle subtended at the body by the X * earth's equator. X */ Xta_par (tha, tdec, phi, ht, ehp, aha, adec) Xdouble tha, tdec, phi, ht, ehp; Xdouble *aha, *adec; X{ X static double last_phi, last_ht, rsp, rcp; X double rp; /* distance to object in Earth radii */ X double ctha; X double stdec, ctdec; X double tdtha, dtha; X double caha; X X /* avoid calcs involving the same phi and ht */ X if (phi != last_phi || ht != last_ht) { X double cphi, sphi, u; X cphi = cos(phi); X sphi = sin(phi); X u = atan(9.96647e-1*sphi/cphi); X rsp = (9.96647e-1*sin(u))+(ht*sphi); X rcp = cos(u)+(ht*cphi); X last_phi = phi; X last_ht = ht; X } X X rp = 1/sin(ehp); X X ctha = cos(tha); X stdec = sin(tdec); X ctdec = cos(tdec); X tdtha = (rcp*sin(tha))/((rp*ctdec)-(rcp*ctha)); X dtha = atan(tdtha); X *aha = tha+dtha; X caha = cos(*aha); X range (aha, 2*PI); X *adec = atan(caha*(rp*stdec-rsp)/(rp*ctdec*ctha-rcp)); X} X X/* given the apparent ha and dec, aha and adec, the geographical latitude, phi, X * the height above sea-level (as a fraction of the earths radius, 6378.16km), X * ht, and the equatorial horizontal parallax, ehp, find the true ha and dec, X * tha and tdec allowing for parallax. X * all angles in radians. ehp is the angle subtended at the body by the X * earth's equator. X * uses ta_par() iteratively: find a set of true ha/dec that converts back X * to the given apparent ha/dec. X */ Xat_par (aha, adec, phi, ht, ehp, tha, tdec) Xdouble aha, adec, phi, ht, ehp; Xdouble *tha, *tdec; X{ X double nha, ndec; /* ha/dec corres. to current true guesses */ X double eha, edec; /* error in ha/dec */ X X /* first guess for true is just the apparent */ X *tha = aha; X *tdec = adec; X X while (1) { X ta_par (*tha, *tdec, phi, ht, ehp, &nha, &ndec); X eha = aha - nha; X edec = adec - ndec; X if (fabs(eha)<1e-6 && fabs(edec)<1e-6) X break; X *tha += eha; X *tdec += edec; X } X} EOFxEOF len=`wc -c < parallax.c` if expr $len != 2280 > /dev/null then echo Length of parallax.c is $len but it should be 2280. fi echo x pelement.c sed -e 's/^X//' << 'EOFxEOF' > pelement.c X#include <stdio.h> X#include <math.h> X#include "astro.h" X X/* this array contains polynomial coefficients to find the various orbital X * elements for the mean orbit at any instant in time for each major planet. X * the first five elements are in the form a0 + a1*t + a2*t**2 + a3*t**3, X * where t is the number of Julian centuries of 36525 Julian days since 1900 X * Jan 0.5. the last three elements are constants. X * X * the orbital element (column) indeces are: X * [ 0- 3]: coefficients for mean longitude, in degrees; X * [ 4- 7]: coefficients for longitude of the perihelion, in degrees; X * [ 8-11]: coefficients for eccentricity; X * [12-15]: coefficients for inclination, in degrees; X * [16-19]: coefficients for longitude of the ascending node, in degrees; X * [20]: semi-major axis, in AU; X * [21]: angular diameter at 1 AU, in arcsec; X * [22]: standard visual magnitude, ie, the visual magnitude of the planet X * when at a distance of 1 AU from both the Sun and the Earth and X * with zero phase angle. X * X * the planent (row) indeces are: X * [0]: Mercury; [1]: Venus; [2]: Mars; [3]: Jupiter; [4]: Saturn; X * [5]: Uranus; [6]: Neptune; [7]: Pluto. X */ X#define NPELE (5*4 + 3) /* 4 coeffs for ea of 5 elems, + 3 constants */ Xstatic double elements[8][NPELE] = { X X { /* mercury... */ X X 178.179078, 415.2057519, 3.011e-4, 0.0, X 75.899697, 1.5554889, 2.947e-4, 0.0, X .20561421, 2.046e-5, 3e-8, 0.0, X 7.002881, 1.8608e-3, -1.83e-5, 0.0, X 47.145944, 1.1852083, 1.739e-4, 0.0, X .3870986, 6.74, -0.42 X }, X X { /* venus... */ X X 342.767053, 162.5533664, 3.097e-4, 0.0, X 130.163833, 1.4080361, -9.764e-4, 0.0, X 6.82069e-3, -4.774e-5, 9.1e-8, 0.0, X 3.393631, 1.0058e-3, -1e-6, 0.0, X 75.779647, .89985, 4.1e-4, 0.0, X .7233316, 16.92, -4.4 X }, X X { /* mars... */ X X 293.737334, 53.17137642, 3.107e-4, 0.0, X 3.34218203e2, 1.8407584, 1.299e-4, -1.19e-6, X 9.33129e-2, 9.2064e-5, 7.7e-8, 0.0, X 1.850333, -6.75e-4, 1.26e-5, 0.0, X 48.786442, .7709917, -1.4e-6, -5.33e-6, X 1.5236883, 9.36, -1.52 X }, X X { /* jupiter... */ X X 238.049257, 8.434172183, 3.347e-4, -1.65e-6, X 1.2720972e1, 1.6099617, 1.05627e-3, -3.43e-6, X 4.833475e-2, 1.6418e-4, -4.676e-7, -1.7e-9, X 1.308736, -5.6961e-3, 3.9e-6, 0.0, X 99.443414, 1.01053, 3.5222e-4, -8.51e-6, X 5.202561, 196.74, -9.4 X }, X X { /* saturn... */ X X 266.564377, 3.398638567, 3.245e-4, -5.8e-6, X 9.1098214e1, 1.9584158, 8.2636e-4, 4.61e-6, X 5.589232e-2, -3.455e-4, -7.28e-7, 7.4e-10, X 2.492519, -3.9189e-3, -1.549e-5, 4e-8, X 112.790414, .8731951, -1.5218e-4, -5.31e-6, X 9.554747, 165.6, -8.88 X }, X X { /* uranus... */ X X 244.19747, 1.194065406, 3.16e-4, -6e-7, X 1.71548692e2, 1.4844328, 2.372e-4, -6.1e-7, X 4.63444e-2, -2.658e-5, 7.7e-8, 0.0, X .772464, 6.253e-4, 3.95e-5, 0.0, X 73.477111, .4986678, 1.3117e-3, 0.0, X 19.21814, 65.8, -7.19 X }, X X { /* neptune... */ X X 84.457994, .6107942056, 3.205e-4, -6e-7, X 4.6727364e1, 1.4245744, 3.9082e-4, -6.05e-7, X 8.99704e-3, 6.33e-6, -2e-9, 0.0, X 1.779242, -9.5436e-3, -9.1e-6, 0.0, X 130.681389, 1.098935, 2.4987e-4, -4.718e-6, X 30.10957, 62.2, -6.87 X }, X X { /* pluto...(osculating 1984 jan 21) */ X X 95.3113544, .3980332167, 0.0, 0.0, X 224.017, 0.0, 0.0, 0.0, X .25515, 0.0, 0.0, 0.0, X 17.1329, 0.0, 0.0, 0.0, X 110.191, 0.0, 0.0, 0.0, X 39.8151, 8.2, -1.0 X } X}; X X/* given a modified Julian date, mjd, return the elements for the mean orbit X * at that instant of all the major planets, together with their X * mean daily motions in longitude, angular diameter and standard visual X * magnitude. X * plan[i][j] contains all the values for all the planets at mjd, such that X * i = 0..7: mercury, venus, mars, jupiter, saturn, unranus, neptune, pluto; X * j = 0..8: mean longitude, mean daily motion in longitude, longitude of X * the perihelion, eccentricity, inclination, longitude of the ascending X * node, length of the semi-major axis, angular diameter from 1 AU, and X * the standard visual magnitude (see elements[][] comment, above). X */ Xpelement (mjd, plan) Xdouble mjd; Xdouble plan[8][9]; X{ X register double *ep, *pp; X register double t = mjd/36525.; X double aa; X int planet, i; X X for (planet = 0; planet < 8; planet++) { X ep = elements[planet]; X pp = plan[planet]; X aa = ep[1]*t; X pp[0] = ep[0] + 360.*(aa-(long)aa) + (ep[3]*t + ep[2])*t*t; X range (pp, 360.); X pp[1] = (ep[1]*9.856263e-3) + (ep[2] + ep[3])/36525; X X for (i = 4; i < 20; i += 4) X pp[i/4+1] = ((ep[i+3]*t + ep[i+2])*t + ep[i+1])*t + ep[i+0]; X X pp[6] = ep[20]; X pp[7] = ep[21]; X pp[8] = ep[22]; X } X} EOFxEOF len=`wc -c < pelement.c` if expr $len != 4797 > /dev/null then echo Length of pelement.c is $len but it should be 4797. fi echo x plans.c sed -e 's/^X//' << 'EOFxEOF' > plans.c X#include <stdio.h> X#include <math.h> X#include "astro.h" X X#define TWOPI (2*PI) X#define mod2PI(x) ((x) - (long)((x)/TWOPI)*TWOPI) X X/* given a modified Julian date, mjd, and a planet, p, find: X * lpd0: heliocentric longitude, X * psi0: heliocentric latitude, X * rp0: distance from the sun to the planet, X * rho0: distance from the Earth to the planet, X * none corrected for light time, ie, they are the true values for the X * given instant. X * lam: geocentric ecliptic longitude, X * bet: geocentric ecliptic latitude, X * each corrected for light time, ie, they are the apparent values as X * seen from the center of the Earth for the given instant. X * dia: angular diameter in arcsec at 1 AU, X * mag: visual magnitude when 1 AU from sun and earth at 0 phase angle. X * X * all angles are in radians, all distances in AU. X * the mean orbital elements are found by calling pelement(), then mutual X * perturbation corrections are applied as necessary. X * X * corrections for nutation and abberation must be made by the caller. The RA X * and DEC calculated from the fully-corrected ecliptic coordinates are then X * the apparent geocentric coordinates. Further corrections can be made, if X * required, for atmospheric refraction and geocentric parallax although the X * intrinsic error herein of about 10 arcseconds is usually the dominant X * error at this stage. X * TODO: combine the several intermediate expressions when get a good compiler. X */ Xplans (mjd, p, lpd0, psi0, rp0, rho0, lam, bet, dia, mag) Xdouble mjd; Xint p; Xdouble *lpd0, *psi0, *rp0, *rho0, *lam, *bet, *dia, *mag; X{ X static double plan[8][9]; X static double lastmjd = -10000; X double dl; /* perturbation correction for longitude */ X double dr; /* " orbital radius */ X double dml; /* " mean longitude */ X double ds; /* " eccentricity */ X double dm; /* " mean anomaly */ X double da; /* " semi-major axis */ X double dhl; /* " heliocentric longitude */ X double lsn, rsn; /* true geocentric longitude of sun and sun-earth rad */ X double mas; /* mean anomaly of the sun */ X double re; /* radius of earth's orbit */ X double lg; /* longitude of earth */ X double map[8]; /* array of mean anomalies for each planet */ X double lpd, psi, rp, rho; X double ll, sll, cll; X double t; X double dt; X int pass; X int j; X double s, ma; X double nu, ea; X double lp, om; X double lo, slo, clo; X double inc, y; X double spsi, cpsi; X double rpd; X X /* only need to fill in plan[] once for a given mjd */ X if (mjd != lastmjd) { X pelement (mjd, plan); X lastmjd = mjd; X } X X dt = 0; X t = mjd/36525.; X sunpos (mjd, &lsn, &rsn); X masun (mjd, &mas); X re = rsn; X lg = lsn+PI; X X /* first find the true position of the planet at mjd. X * then repeat a second time for a slightly different time based X * on the position found in the first pass to account for light-travel X * time. X */ X for (pass = 0; pass < 2; pass++) { X X for (j = 0; j < 8; j++) X map[j] = degrad(plan[j][0]-plan[j][2]-dt*plan[j][1]); X X /* set initial corrections to 0. X * then modify as necessary for the planet of interest. X */ X dl = 0; X dr = 0; X dml = 0; X ds = 0; X dm = 0; X da = 0; X dhl = 0; X X switch (p) { X X case MERCURY: X p_mercury (map, &dl, &dr); X break; X X case VENUS: X p_venus (t, mas, map, &dl, &dr, &dml, &dm); X break; X X case MARS: X p_mars (mas, map, &dl, &dr, &dml, &dm); X break; X X case JUPITER: X p_jupiter (t, plan[p][3], &dml, &ds, &dm, &da); X break; X X case SATURN: X p_saturn (t, plan[p][3], &dml, &ds, &dm, &da, &dhl); X break; X X case URANUS: X p_uranus (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl); X break; X X case NEPTUNE: X p_neptune (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl); X break; X X case PLUTO: X /* no perturbation theory for pluto */ X break; X } X X s = plan[p][3]+ds; X ma = map[p]+dm; X anomaly (ma, s, &nu, &ea); X rp = (plan[p][6]+da)*(1-s*s)/(1+s*cos(nu)); X lp = raddeg(nu)+plan[p][2]+raddeg(dml-dm); X lp = degrad(lp); X om = degrad(plan[p][5]); X lo = lp-om; X slo = sin(lo); X clo = cos(lo); X inc = degrad(plan[p][4]); X rp = rp+dr; X spsi = slo*sin(inc); X y = slo*cos(inc); X psi = asin(spsi)+dhl; X spsi = sin(psi); X lpd = atan(y/clo)+om+degrad(dl); X if (clo<0) lpd += PI; X range (&lpd, TWOPI); X cpsi = cos(psi); X rpd = rp*cpsi; X ll = lpd-lg; X rho = sqrt(re*re+rp*rp-2*re*rp*cpsi*cos(ll)); X X /* when we view a planet we see it in the position it occupied X * dt days ago, where rho is the distance between it and earth, X * in AU. use this as the new time for the next pass. X */ X dt = rho*5.775518e-3; X X if (pass == 0) { X /* save heliocentric coordinates after first pass since, being X * true, they are NOT to be corrected for light-travel time. X */ X *lpd0 = lpd; X range (lpd0, TWOPI); X *psi0 = psi; X *rp0 = rp; X *rho0 = rho; X } X } X X sll = sin(ll); X cll = cos(ll); X if (p < MARS) X *lam = atan(-1*rpd*sll/(re-rpd*cll))+lg+PI; X else X *lam = atan(re*sll/(rpd-re*cll))+lpd; X range (lam, TWOPI); X *bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*re*sll)); X *dia = plan[p][7]; X *mag = plan[p][8]; X} X X/* set auxilliary variables used for jupiter, saturn, uranus, and neptune */ Xstatic Xaux_jsun (t, x1, x2, x3, x4, x5, x6) Xdouble t; Xdouble *x1, *x2, *x3, *x4, *x5, *x6; X{ X *x1 = t/5+0.1; X *x2 = mod2PI(4.14473+5.29691e1*t); X *x3 = mod2PI(4.641118+2.132991e1*t); X *x4 = mod2PI(4.250177+7.478172*t); X *x5 = 5 * *x3 - 2 * *x2; X *x6 = 2 * *x2 - 6 * *x3 + 3 * *x4; X} X X/* find the mean anomaly of the sun at mjd. X * this is the same as that used in sun() but when it was converted to C it X * was not known it would be required outside that routine. X * TODO: add an argument to sun() to return mas and eliminate this routine. X */ Xstatic Xmasun (mjd, mas) Xdouble mjd; Xdouble *mas; X{ X double t, t2; X double a, b; X X t = mjd/36525; X t2 = t*t; X a = 9.999736042e1*t; X b = 360.*(a-(long)a); X *mas = degrad (3.5847583e2-(1.5e-4+3.3e-6*t)*t2+b); X} X X/* perturbations for mercury */ Xstatic Xp_mercury (map, dl, dr) Xdouble map[]; Xdouble *dl, *dr; X{ X *dl = 2.04e-3*cos(5*map[2-1]-2*map[1-1]+2.1328e-1)+ X 1.03e-3*cos(2*map[2-1]-map[1-1]-2.8046)+ X 9.1e-4*cos(2*map[3]-map[1-1]-6.4582e-1)+ X 7.8e-4*cos(5*map[2-1]-3*map[1-1]+1.7692e-1); X X *dr = 7.525e-6*cos(2*map[3]-map[1-1]+9.25251e-1)+ X 6.802e-6*cos(5*map[2-1]-3*map[1-1]-4.53642)+ X 5.457e-6*cos(2*map[2-1]-2*map[1-1]-1.24246)+ X 3.569e-6*cos(5*map[2-1]-map[1-1]-1.35699); X} X X/* ....venus */ Xstatic Xp_venus (t, mas, map, dl, dr, dml, dm) Xdouble t, mas, map[]; Xdouble *dl, *dr, *dml, *dm; X{ X *dml = degrad (7.7e-4*sin(4.1406+t*2.6227)); X *dm = *dml; X X *dl = 3.13e-3*cos(2*mas-2*map[2-1]-2.587)+ X 1.98e-3*cos(3*mas-3*map[2-1]+4.4768e-2)+ X 1.36e-3*cos(mas-map[2-1]-2.0788)+ X 9.6e-4*cos(3*mas-2*map[2-1]-2.3721)+ X 8.2e-4*cos(map[3]-map[2-1]-3.6318); X X *dr = 2.2501e-5*cos(2*mas-2*map[2-1]-1.01592)+ X 1.9045e-5*cos(3*mas-3*map[2-1]+1.61577)+ X 6.887e-6*cos(map[3]-map[2-1]-2.06106)+ X 5.172e-6*cos(mas-map[2-1]-5.08065e-1)+ X 3.62e-6*cos(5*mas-4*map[2-1]-1.81877)+ X 3.283e-6*cos(4*mas-4*map[2-1]+1.10851)+ X 3.074e-6*cos(2*map[3]-2*map[2-1]-9.62846e-1); X} X X/* ....mars */ Xstatic Xp_mars (mas, map, dl, dr, dml, dm) Xdouble mas, map[]; Xdouble *dl, *dr, *dml, *dm; X{ X double a; X X a = 3*map[3]-8*map[2]+4*mas; X *dml = degrad (-1*(1.133e-2*sin(a)+9.33e-3*cos(a))); X *dm = *dml; X X *dl = 7.05e-3*cos(map[3]-map[2]-8.5448e-1)+ X 6.07e-3*cos(2*map[3]-map[2]-3.2873)+ X 4.45e-3*cos(2*map[3]-2*map[2]-3.3492)+ X 3.88e-3*cos(mas-2*map[2]+3.5771e-1)+ X 2.38e-3*cos(mas-map[2]+6.1256e-1)+ X 2.04e-3*cos(2*mas-3*map[2]+2.7688)+ X 1.77e-3*cos(3*map[2]-map[2-1]-1.0053)+ X 1.36e-3*cos(2*mas-4*map[2]+2.6894)+ X 1.04e-3*cos(map[3]+3.0749e-1); X X *dr = 5.3227e-5*cos(map[3]-map[2]+7.17864e-1)+ X 5.0989e-5*cos(2*map[3]-2*map[2]-1.77997)+ X 3.8278e-5*cos(2*map[3]-map[2]-1.71617)+ X 1.5996e-5*cos(mas-map[2]-9.69618e-1)+ X 1.4764e-5*cos(2*mas-3*map[2]+1.19768)+ X 8.966e-6*cos(map[3]-2*map[2]+7.61225e-1); X *dr += 7.914e-6*cos(3*map[3]-2*map[2]-2.43887)+ X 7.004e-6*cos(2*map[3]-3*map[2]-1.79573)+ X 6.62e-6*cos(mas-2*map[2]+1.97575)+ X 4.93e-6*cos(3*map[3]-3*map[2]-1.33069)+ X 4.693e-6*cos(3*mas-5*map[2]+3.32665)+ X 4.571e-6*cos(2*mas-4*map[2]+4.27086)+ X 4.409e-6*cos(3*map[3]-map[2]-2.02158); X} X X/* ....jupiter */ Xstatic Xp_jupiter (t, s, dml, ds, dm, da) Xdouble t, s; Xdouble *dml, *ds, *dm, *da; X{ X double dp; X double x1, x2, x3, x4, x5, x6, x7; X double sx3, cx3, s2x3, c2x3; X double sx5, cx5, s2x5; X double sx6; X double sx7, cx7, s2x7, c2x7, s3x7, c3x7, s4x7, c4x7, c5x7; X X aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6); X x7 = x3-x2; X sx3 = sin(x3); X cx3 = cos(x3); X s2x3 = sin(2*x3); X c2x3 = cos(2*x3); X sx5 = sin(x5); X cx5 = cos(x5); X s2x5 = sin(2*x5); X sx6 = sin(x6); X sx7 = sin(x7); X cx7 = cos(x7); X s2x7 = sin(2*x7); X c2x7 = cos(2*x7); X s3x7 = sin(3*x7); X c3x7 = cos(3*x7); X s4x7 = sin(4*x7); X c4x7 = cos(4*x7); X c5x7 = cos(5*x7); X X *dml = (3.31364e-1-(1.0281e-2+4.692e-3*x1)*x1)*sx5+ X (3.228e-3-(6.4436e-2-2.075e-3*x1)*x1)*cx5- X (3.083e-3+(2.75e-4-4.89e-4*x1)*x1)*s2x5+ X 2.472e-3*sx6+1.3619e-2*sx7+1.8472e-2*s2x7+6.717e-3*s3x7+ X 2.775e-3*s4x7+6.417e-3*s2x7*sx3+ X (7.275e-3-1.253e-3*x1)*sx7*sx3+ X 2.439e-3*s3x7*sx3-(3.5681e-2+1.208e-3*x1)*sx7*cx3; X *dml += -3.767e-3*c2x7*sx3-(3.3839e-2+1.125e-3*x1)*cx7*sx3- X 4.261e-3*s2x7*cx3+ X (1.161e-3*x1-6.333e-3)*cx7*cx3+ X 2.178e-3*cx3-6.675e-3*c2x7*cx3-2.664e-3*c3x7*cx3- X 2.572e-3*sx7*s2x3-3.567e-3*s2x7*s2x3+2.094e-3*cx7*c2x3+ X 3.342e-3*c2x7*c2x3; X *dml = degrad(*dml); X X *ds = (3606+(130-43*x1)*x1)*sx5+(1289-580*x1)*cx5-6764*sx7*sx3- X 1110*s2x7*sx3-224*s3x7*sx3-204*sx3+(1284+116*x1)*cx7*sx3+ X 188*c2x7*sx3+(1460+130*x1)*sx7*cx3+224*s2x7*cx3-817*cx3+ X 6074*cx3*cx7+992*c2x7*cx3+ X 508*c3x7*cx3+230*c4x7*cx3+108*c5x7*cx3; X *ds += -(956+73*x1)*sx7*s2x3+448*s2x7*s2x3+137*s3x7*s2x3+ X (108*x1-997)*cx7*s2x3+480*c2x7*s2x3+148*c3x7*s2x3+ X (99*x1-956)*sx7*c2x3+490*s2x7*c2x3+ X 158*s3x7*c2x3+179*c2x3+(1024+75*x1)*cx7*c2x3- X 437*c2x7*c2x3-132*c3x7*c2x3; X *ds *= 1e-7; X X dp = (7.192e-3-3.147e-3*x1)*sx5-4.344e-3*sx3+ X (x1*(1.97e-4*x1-6.75e-4)-2.0428e-2)*cx5+ X 3.4036e-2*cx7*sx3+(7.269e-3+6.72e-4*x1)*sx7*sx3+ X 5.614e-3*c2x7*sx3+2.964e-3*c3x7*sx3+3.7761e-2*sx7*cx3+ X 6.158e-3*s2x7*cx3- X 6.603e-3*cx7*cx3-5.356e-3*sx7*s2x3+2.722e-3*s2x7*s2x3+ X 4.483e-3*cx7*s2x3-2.642e-3*c2x7*s2x3+4.403e-3*sx7*c2x3- X 2.536e-3*s2x7*c2x3+5.547e-3*cx7*c2x3-2.689e-3*c2x7*c2x3; X X *dm = *dml-(degrad(dp)/s); X X *da = 205*cx7-263*cx5+693*c2x7+312*c3x7+147*c4x7+299*sx7*sx3+ X 181*c2x7*sx3+204*s2x7*cx3+111*s3x7*cx3-337*cx7*cx3- X 111*c2x7*cx3; X *da *= 1e-6; X} X X/* ....saturn */ Xstatic Xp_saturn (t, s, dml, ds, dm, da, dhl) Xdouble t, s; Xdouble *dml, *ds, *dm, *da, *dhl; X{ X double dp; X double x1, x2, x3, x4, x5, x6, x7, x8; X double sx3, cx3, s2x3, c2x3, s3x3, c3x3, s4x3, c4x3; X double sx5, cx5, s2x5, c2x5; X double sx6; X double sx7, cx7, s2x7, c2x7, s3x7, c3x7, s4x7, c4x7, c5x7, s5x7; X double s2x8, c2x8, s3x8, c3x8; X X aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6); X x7 = x3-x2; X sx3 = sin(x3); X cx3 = cos(x3); X s2x3 = sin(2*x3); X c2x3 = cos(2*x3); X sx5 = sin(x5); X cx5 = cos(x5); X s2x5 = sin(2*x5); X sx6 = sin(x6); X sx7 = sin(x7); X cx7 = cos(x7); X s2x7 = sin(2*x7); X c2x7 = cos(2*x7); X s3x7 = sin(3*x7); X c3x7 = cos(3*x7); X s4x7 = sin(4*x7); X c4x7 = cos(4*x7); X c5x7 = cos(5*x7); X X s3x3 = sin(3*x3); X c3x3 = cos(3*x3); X s4x3 = sin(4*x3); X c4x3 = cos(4*x3); X c2x5 = cos(2*x5); X s5x7 = sin(5*x7); X x8 = x4-x3; X s2x8 = sin(2*x8); X c2x8 = cos(2*x8); X s3x8 = sin(3*x8); X c3x8 = cos(3*x8); X X *dml = 7.581e-3*s2x5-7.986e-3*sx6-1.48811e-1*sx7-4.0786e-2*s2x7- X (8.14181e-1-(1.815e-2-1.6714e-2*x1)*x1)*sx5- X (1.0497e-2-(1.60906e-1-4.1e-3*x1)*x1)*cx5-1.5208e-2*s3x7- X 6.339e-3*s4x7-6.244e-3*sx3-1.65e-2*s2x7*sx3+ X (8.931e-3+2.728e-3*x1)*sx7*sx3-5.775e-3*s3x7*sx3+ X (8.1344e-2+3.206e-3*x1)*cx7*sx3+1.5019e-2*c2x7*sx3; X *dml += (8.5581e-2+2.494e-3*x1)*sx7*cx3+1.4394e-2*c2x7*cx3+ X (2.5328e-2-3.117e-3*x1)*cx7*cx3+ X 6.319e-3*c3x7*cx3+6.369e-3*sx7*s2x3+9.156e-3*s2x7*s2x3+ X 7.525e-3*s3x8*s2x3-5.236e-3*cx7*c2x3-7.736e-3*c2x7*c2x3- X 7.528e-3*c3x8*c2x3; X *dml = degrad(*dml); X X *ds = (-7927+(2548+91*x1)*x1)*sx5+(13381+(1226-253*x1)*x1)*cx5+ X (248-121*x1)*s2x5-(305+91*x1)*c2x5+412*s2x7+12415*sx3+ X (390-617*x1)*sx7*sx3+(165-204*x1)*s2x7*sx3+26599*cx7*sx3- X 4687*c2x7*sx3-1870*c3x7*sx3-821*c4x7*sx3- X 377*c5x7*sx3+497*c2x8*sx3+(163-611*x1)*cx3; X *ds += -12696*sx7*cx3-4200*s2x7*cx3-1503*s3x7*cx3-619*s4x7*cx3- X 268*s5x7*cx3-(282+1306*x1)*cx7*cx3+(-86+230*x1)*c2x7*cx3+ X 461*s2x8*cx3-350*s2x3+(2211-286*x1)*sx7*s2x3- X 2208*s2x7*s2x3-568*s3x7*s2x3-346*s4x7*s2x3- X (2780+222*x1)*cx7*s2x3+(2022+263*x1)*c2x7*s2x3+248*c3x7*s2x3+ X 242*s3x8*s2x3+467*c3x8*s2x3-490*c2x3-(2842+279*x1)*sx7*c2x3; X *ds += (128+226*x1)*s2x7*c2x3+224*s3x7*c2x3+ X (-1594+282*x1)*cx7*c2x3+(2162-207*x1)*c2x7*c2x3+ X 561*c3x7*c2x3+343*c4x7*c2x3+469*s3x8*c2x3-242*c3x8*c2x3- X 205*sx7*s3x3+262*s3x7*s3x3+208*cx7*c3x3-271*c3x7*c3x3- X 382*c3x7*s4x3-376*s3x7*c4x3; X *ds *= 1e-7; X X dp = (7.7108e-2+(7.186e-3-1.533e-3*x1)*x1)*sx5-7.075e-3*sx7+ X (4.5803e-2-(1.4766e-2+5.36e-4*x1)*x1)*cx5-7.2586e-2*cx3- X 7.5825e-2*sx7*sx3-2.4839e-2*s2x7*sx3-8.631e-3*s3x7*sx3- X 1.50383e-1*cx7*cx3+2.6897e-2*c2x7*cx3+1.0053e-2*c3x7*cx3- X (1.3597e-2+1.719e-3*x1)*sx7*s2x3+1.1981e-2*s2x7*c2x3; X dp += -(7.742e-3-1.517e-3*x1)*cx7*s2x3+ X (1.3586e-2-1.375e-3*x1)*c2x7*c2x3- X (1.3667e-2-1.239e-3*x1)*sx7*c2x3+ X (1.4861e-2+1.136e-3*x1)*cx7*c2x3- X (1.3064e-2+1.628e-3*x1)*c2x7*c2x3; X X *dm = *dml-(degrad(dp)/s); X X *da = 572*sx5-1590*s2x7*cx3+2933*cx5-647*s3x7*cx3+33629*cx7- X 344*s4x7*cx3-3081*c2x7+2885*cx7*cx3-1423*c3x7+ X (2172+102*x1)*c2x7*cx3-671*c4x7+296*c3x7*cx3-320*c5x7- X 267*s2x7*s2x3+1098*sx3-778*cx7*s2x3-2812*sx7*sx3; X *da += 495*c2x7*s2x3+688*s2x7*sx3+250*c3x7*s2x3-393*s3x7*sx3- X 856*sx7*c2x3-228*s4x7*sx3+441*s2x7*c2x3+2138*cx7*sx3+ X 296*c2x7*c2x3-999*c2x7*sx3+211*c3x7*c2x3-642*c3x7*sx3- X 427*sx7*s3x3-325*c4x7*sx3+398*s3x7*s3x3-890*cx3+ X 344*cx7*c3x3+2206*sx7*cx3-427*c3x7*c3x3; X *da *= 1e-6; X X *dhl = 7.47e-4*cx7*sx3+1.069e-3*cx7*cx3+2.108e-3*s2x7*s2x3+ X 1.261e-3*c2x7*s2x3+1.236e-3*s2x7*c2x3-2.075e-3*c2x7*c2x3; X *dhl = degrad(*dhl); X} X X/* ....uranus */ Xstatic Xp_uranus (t, s, dl, dr, dml, ds, dm, da, dhl) Xdouble t, s; Xdouble *dl, *dr, *dml, *ds, *dm, *da, *dhl; X{ X double dp; X double x1, x2, x3, x4, x5, x6; X double x8, x9, x10, x11, x12; X double sx4, cx4, s2x4, c2x4; X double sx9, cx9, s2x9, c2x9; X double sx11, cx11; X X aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6); X X x8 = mod2PI(1.46205+3.81337*t); X x9 = 2*x8-x4; X sx9 = sin(x9); X cx9 = cos(x9); X s2x9 = sin(2*x9); X c2x9 = cos(2*x9); X X x10 = x4-x2; X x11 = x4-x3; X x12 = x8-x4; X X *dml = (8.64319e-1-1.583e-3*x1)*sx9+(8.2222e-2-6.833e-3*x1)*cx9+ X 3.6017e-2*s2x9-3.019e-3*c2x9+8.122e-3*sin(x6); X *dml = degrad(*dml); X X dp = 1.20303e-1*sx9+6.197e-3*s2x9+(1.9472e-2-9.47e-4*x1)*cx9; X *dm = *dml-(degrad(dp)/s); X X *ds = (163*x1-3349)*sx9+20981*cx9+1311*c2x9; X *ds *= 1e-7; X X *da = -3.825e-3*cx9; X X *dl = (1.0122e-2-9.88e-4*x1)*sin(x4+x11)+ X (-3.8581e-2+(2.031e-3-1.91e-3*x1)*x1)*cos(x4+x11)+ X (3.4964e-2-(1.038e-3-8.68e-4*x1)*x1)*cos(2*x4+x11)+ X 5.594e-3*sin(x4+3*x12)-1.4808e-2*sin(x10)- X 5.794e-3*sin(x11)+2.347e-3*cos(x11)+9.872e-3*sin(x12)+ X 8.803e-3*sin(2*x12)-4.308e-3*sin(3*x12); X X sx11 = sin(x11); X cx11 = cos(x11); X sx4 = sin(x4); X cx4 = cos(x4); X s2x4 = sin(2*x4); X c2x4 = cos(2*x4); X *dhl = (4.58e-4*sx11-6.42e-4*cx11-5.17e-4*cos(4*x12))*sx4- X (3.47e-4*sx11+8.53e-4*cx11+5.17e-4*sin(4*x11))*cx4+ X 4.03e-4*(cos(2*x12)*s2x4+sin(2*x12)*c2x4); X *dhl = degrad(*dhl); X X *dr = -25948+4985*cos(x10)-1230*cx4+3354*cos(x11)+904*cos(2*x12)+ X 894*(cos(x12)-cos(3*x12))+(5795*cx4-1165*sx4+1388*c2x4)*sx11+ X (1351*cx4+5702*sx4+1388*s2x4)*cos(x11); X *dr *= 1e-6; X} X X/* ....neptune */ Xstatic Xp_neptune (t, s, dl, dr, dml, ds, dm, da, dhl) Xdouble t, s; Xdouble *dl, *dr, *dml, *ds, *dm, *da, *dhl; X{ X double dp; X double x1, x2, x3, x4, x5, x6; X double x8, x9, x10, x11, x12; X double sx8, cx8; X double sx9, cx9, s2x9, c2x9; X double s2x12, c2x12; X X aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6); X X x8 = mod2PI(1.46205+3.81337*t); X x9 = 2*x8-x4; X sx9 = sin(x9); X cx9 = cos(x9); X s2x9 = sin(2*x9); X c2x9 = cos(2*x9); X X x10 = x8-x2; X x11 = x8-x3; X x12 = x8-x4; X X *dml = (1.089e-3*x1-5.89833e-1)*sx9+(4.658e-3*x1-5.6094e-2)*cx9- X 2.4286e-2*s2x9; X *dml = degrad(*dml); X X dp = 2.4039e-2*sx9-2.5303e-2*cx9+6.206e-3*s2x9-5.992e-3*c2x9; X X *dm = *dml-(degrad(dp)/s); X X *ds = 4389*sx9+1129*s2x9+4262*cx9+1089*c2x9; X *ds *= 1e-7; X X *da = 8189*cx9-817*sx9+781*c2x9; X *da *= 1e-6; X X s2x12 = sin(2*x12); X c2x12 = cos(2*x12); X sx8 = sin(x8); X cx8 = cos(x8); X *dl = -9.556e-3*sin(x10)-5.178e-3*sin(x11)+2.572e-3*s2x12- X 2.972e-3*c2x12*sx8-2.833e-3*s2x12*cx8; X X *dhl = 3.36e-4*c2x12*sx8+3.64e-4*s2x12*cx8; X *dhl = degrad(*dhl); X X *dr = -40596+4992*cos(x10)+2744*cos(x11)+2044*cos(x12)+1051*c2x12; X *dr *= 1e-6; X} X EOFxEOF len=`wc -c < plans.c` if expr $len != 17648 > /dev/null then echo Length of plans.c is $len but it should be 17648. fi